Keen - jl

This code simulates the Keen model as specified in Grasselli and Costa Lima (2012)
Authors: B. Costa Lima and M. Grasselli

Links & source

Grasselli, M. R., and B. Costa Lima. 2012. “An Analysis of the Keen Model for Credit Expansion, Asset Price Bubbles and Financial Fragility.” Mathematics and Financial Economics 6 (3): 191–210. https://doi.org/10.1007/s11579-012-0071-8.

https://ms.mcmaster.ca/~grasselli/GrasselliCostaLima_MAFE_online.pdf
https://www.oecd.org/naec/projects/naecinnovationlabevents/matheus_grasselli_OECD_masterclass_2020.pdf
https://ms.mcmaster.ca/~grasselli/keen.html

https://diffeq.sciml.ai/stable/tutorials/ode_example/

• Grasselli, M. R., and B. Costa Lima. 2012. “An Analysis of the Keen Model for Credit Expansion, Asset Price Bubbles and Financial Fragility.” Mathematics and Financial Economics 6 (3): 191–210. https://doi.org/10.1007/s11579-012-0071-8.
• Huu, A. Nguyen, and Bernardo Costa-Lima. 2014. “Orbits in a Stochastic Goodwin–Lotka–Volterra Model.” Journal of Mathematical Analysis and Applications 419 (1): 48–67. http://www.sciencedirect.com/science/article/pii/S0022247X14003862.
• Keen, Steve. 1995. “Finance and Economic Breakdown: Modeling Minsky’s ‘Financial Instability Hypothesis.’” Journal of Post Keynesian Economics 17 (4): 607–35. http://econpapers.repec.org/article/mespostke/v_3a17_3ay_3a1995_3ai_3a4_3ap_3a607-635.htm.
• ———. 2013. “A Monetary Minsky Model of the Great Moderation and the Great Recession.” Journal of Economic Behavior & Organization 86 (February): 221–35. https://doi.org/10.1016/j.jebo.2011.01.010.

Preamble

"C:\Users\ibuckley\OneDrive - CSTO\Code\Grasselli\2012 Keen\"

# using Pkg; 
# Pkg.add("ParameterizedFunctions")
# Pkg.add("DifferentialEquations")
# Pkg.add("PlotlyBase")

using DifferentialEquations # for ODEProblem
using DataFrames
using Plots
plotly() # Choose Plotly as the backend to Plots

┌ Info: Precompiling DataFrames [a93c6f00-e57d-5684-b7b6-d8193f3e46c0]
└ @ Base loading.jl:1260

Plots.PlotlyBackend()

#using Plotly
#Plotly.signin("IanBuckley","Z4HJZpdeW4wss1tnizq8")

Keen model

Notation

• $D$ the amount of debt in real terms
• $r$ is a constant real interest rate and
• $d = D/Y$ is the debt ratio in the economy

Capital stock

• The key insight provided by Minsky is that current cash-flows validate past liabilities and form the basis for future ones.
• In other words, high net profits lead to more borrowing whereas low net profits (possibly negative) lead to a deleveraging of the economy.
• Keen formalizes this insight by taking the change in capital stock to be

$$ \dot{K}=\kappa(1-\omega-r d) Y-\delta K $$

where the rate of new investment is a nonlinear increasing function $\kappa$ of the net profit share $π = (1 − ω − rd)$ and $δ$ is a constant depreciation rate as before.

Total output

$$ \frac{\dot{Y}}{Y}=\frac{\kappa(1-\omega-r d)}{v}-\delta:=g(\omega, d) $$

Employment rate

$$ \frac{\dot{\lambda}}{\lambda}=\frac{\kappa(1-\omega-r d)}{v}-\alpha-\beta-\delta $$

Debt

The new dynamic variable in this model is the amount of debt, which changes based on the difference between new investment and net profits.

$$ \dot{D}=\kappa(1-\omega-r d) Y-(1-\omega-r d) Y $$

$$ \frac{\dot{d}}{d}=\frac{\dot{D}}{D}-\frac{\dot{Y}}{Y}=\frac{\kappa(1-\omega-r d)-(1-\omega-r d)}{d}-\frac{\kappa(1-\omega-r d)}{v}+\delta . $$

Parameters

State variables

• $ω$ - wage share of the economy
• $λ$ - employment rate
• $d$ - debt (or reciprocal debt: $u=1/d$)
• $Y$ - output

Fundamental economic constants Equation (24)

• 𝜈 = 3 capital-to-output ratio
• 𝛼 = 0.025 productivity growth rate
• 𝛽 = 0.02 population growth rate
• 𝛿 = 0.01 depreciation rate

ν = 3            # capital-to-output ratio
α = 0.025        # productivity growth rate
β = 0.02         # population growth rate
δ = 0.01;        # depreciation rate

Phillips curve - Equation (25)

$$\Phi(\lambda)=\frac{\phi_{1}}{(1-\lambda)^{2}}-\phi_{0}$$

ϕ0 = 0.04/(1-0.04^2)        # Phillips curve parameter,  Phi(0.96)=0
ϕ1 = 0.04^3/(1-0.04^2);      # Phillips curve parameter, Phi(0)=-0.04

Net profit share $\pi$

$$\pi = (1-\omega - r d)$$

Investment function - Equation (64)

$$\kappa(x)=\kappa_{0}+\kappa_{1} e^{\kappa_{2} x} $$

$$\kappa_{0}=-0.0065, \quad \kappa_{1}=e^{-5}, \quad \kappa_{2}=20$$

r = 0.03                    # real interest rate  κ
κ0 = -0.0065                # investment function parameter
κ1 = exp(-5)                # investment function parameter
κ2 = 20;                    # investment function parameter

Define Phillips curve $\Phi$, & investment function $\kappa$

Phillips curve - Equation (25)

Φ(λ) = ϕ1./(1-λ).^2- ϕ0              # Phillips curve, equation (25)
Φ′(λ) = 2*ϕ1./(1-λ).^3               # derivative of Phillips curve
Φ_inv(λ) = 1 - (ϕ1/(λ+ϕ0)).^(1/2);    # inverse of Phillips curve

Investment function - Equation (64)

κ(x) = κ0+κ1.*exp(κ2.*x)             # Investment function, equation (64)
κ′(x) = κ1.*κ2.*exp(κ2.*x)           # Derivative of investment function
κ_inv(x) = log((x-κ0)./κ1)./κ2;      # Inverse of investment function

Numerical Values of Interest

• investment at interior equilibrium
• profit at interior equilibrium, equation (41)
• debt ratio at interior equilibirum, equation (42)
• wage share at interior equilibrium, equation (42)
• employment at interior equilibirum, equation (42)
• solution to equation (40) given in equation (66)
• eigenvalues of the Jacobian at the interior, first d0 & second d0 equilibrium

https://www.mathworks.com/help/symbolic/vpasolve.html
Symbolic Math Toolbox

Equation (35)

$$\begin{array}{l} \dot{\omega}=\omega[\Phi(\lambda)-\alpha] \\ \dot{\lambda}=\lambda\left[\frac{\kappa(1-\omega-r d)}{v}-\alpha-\beta-\delta\right] \\ \dot{d}=d\left[r-\frac{\kappa(1-\omega-r d)}{v}+\delta\right]+\kappa(1-\omega-r d)-(1-\omega) \end{array}$$

Equation 40

$$d\left[r-\frac{\kappa(1-r d)}{v}+\delta\right]+\kappa(1-r d)-1=0$$

Equation (42) - equilibrium point

$$ \bar{\pi}_{1}=\kappa^{-1}(\nu(\alpha+\beta+\delta)) $$

$$ \begin{array}{l} \bar{\omega}_{1}=1-\bar{\pi}_{1}-r \frac{v(\alpha+\beta+\delta)-\bar{\pi}_{1}}{\alpha+\beta} \\ \bar{\lambda}_{1}=\Phi^{-1}(\alpha) \\ \bar{d}_{1}=\frac{v(\alpha+\beta+\delta)-\bar{\pi}_{1}}{\alpha+\beta} \end{array} $$

Equation (47)

$$\begin{array}{l} \dot{\omega}=\omega[\Phi(\lambda)-\alpha] \\ \dot{\lambda}=\lambda\left[\frac{\kappa(1-\omega-r/u)}{v}-\alpha-\beta-\delta\right] \\ \dot{d}=d\left[r-\frac{\kappa(1-\omega-r/u)}{v}+\delta\right]+\kappa(1-\omega-r/u)-(1-\omega) \end{array}$$

Interior equilibrium point expressions

See Equation (42)

κ_eq = ν*(α+β+δ)                 # investment at interior equilibrium 
π_eq = κ_inv(ν*(α+β+δ))            # profit at interior equilibrium, equation (41)
d_eq = (κ(π_eq)-π_eq)/(α+β)      # debt ratio at interior equilibirum, equation (42)
ω_eq = 1-π_eq-r*d_eq             # wage share at interior equilibrium, equation (42)
λ_eq = Φ_inv(α);                    # employment at interior equilibirum, equation (42)

# TO DO!
#d_sols(1) = vpasolve(d*(r+delta-fun_kappa(1-r*d)/nu)+fun_kappa(1-r*d)-1,d) ;        # solution to equation (40) given in equation (66)
#d_sols(2) = vpasolve(d*(r+delta-fun_kappa(1-r*d)/nu)+fun_kappa(1-r*d)-1,[20 40]) ;  # solution to equation (40) given in equation (66)
#d_sols

Equation (59)

$r\left[\frac{\kappa^{\prime}\left(\bar{\pi}_{1}\right)}{v}\left(\bar{\pi}_{1}-\kappa\left(\bar{\pi}_{1}\right)+v(\alpha+\beta)\right)-(\alpha+\beta)\right]>0$

Jacobian ~ equation (59)

$$ J(\omega, \lambda, u)=\left[\begin{array}{ccc} \Phi(\lambda)-\alpha & \omega \Phi^{\prime}(\lambda) & 0 \\ -\frac{\lambda \kappa^{\prime}(\pi)}{\nu} & \frac{\kappa(\pi)-v(\alpha+\beta+\delta)}{v} & -\frac{r \lambda \kappa^{\prime}(\pi)}{u^{2} v} \\ \frac{\left(v u^{2}-u\right) \kappa^{\prime}(\pi)-v u^{2}}{v} & 0 & \frac{\kappa(\pi)(1-2 u)+r \kappa^{\prime}(\pi)(1 / u-1)+2 u v(1-\omega)-v(r+\delta)}{v} \end{array}\right] $$

J(𝜔,𝜆,𝑢)= [Φ(𝜆)-𝛼 𝜔.*Φ′(𝜆) 0;
    -𝜆.*𝜅′(1-𝜔-r*z)/𝜈 (𝜅(1-𝜔-r*z)-𝜈*(𝛼+𝛽+𝛿))/𝜈 -r*𝜆.*𝜅′(1-𝜔-r*z)/𝜈;
    ((z-𝜈).*𝜅′(1-𝜔-r*z)+𝜈)/𝜈 0 (𝜈*(r+𝛿)-𝜅(1-𝜔-r*z)+r*(z-𝜈).*𝜅′(1-𝜔-r*z))/𝜈];

# TO DO!!!
#E1=eig(J(omega_eq,lambda_eq,d_eq)) # eigenvalues of the Jacobian at the interior equilibrium
#E2=eig(J(0,0,d_sols(1)))           # eigenvalues of the Jacobian at the first d0 equilibrium
#E3=eig(J(0,0,d_sols(2)))           # eigenvalues of the Jacobian at the second d0 equilibrium

cond59= κ′(π_eq)/ν*(π_eq-ν*δ)-(α+β) # expression in brackets in condition (59)

0.10573866662583055

Sample paths of the Keen model using DifferentialEquations

https://diffeq.sciml.ai/stable/tutorials/ode_example/

function keen!(dz,z,p,t) 
    log_ω, tan_λ, π_n = z                         # state variables
    α, β, δ, ν, r, (ϕ0,ϕ1), (κ0,κ1,κ2) = p         # parameters & functions κ, Φ ??? ϕ0,ϕ1,κ0,κ1,κ2
                                                  # Recover λ & ω
    λ = atan(tan_λ)/π + 0.5
    ω = exp(log_ω)
    
    # Functions 
    Φ(λ) = ϕ1./(1-λ).^2 - ϕ0              # Phillips curve, equation (25)
    κ(x) = κ0 + κ1.*exp(κ2.*x)             # Investment function, equation (64)
    
    g_Y = κ(π_n)/ν - δ           # Output function (31)
   
    # Equation (47)
    dz[1] = dlog_ω = Φ(λ) - α                                         #d(log_ω)/dt
    dz[2] = dtan_λ = (1 + tan_λ^2)*π*λ*(g_Y - α - β)                  #d(tan_λ)/dt (constant π!!)
    dz[3] = dπ_n = -ω*dlog_ω - r*(κ(π_n) - π_n) + (1 - ω - π_n)*g_Y   #d(π_n)/dt
end

keen! (generic function with 1 method)

convert

https://ms.mcmaster.ca/~grasselli/keen.html

function convert(state_variables,r)
    n = length(state_variables)
    if n==2
        ω, λ = state_variables
    elseif n==3
        ω, λ, d = state_variables
    else
        ω, λ, d, p = state_variables
    end
    
    log_ω = log.(ω)
    tan_λ = tan.((λ-0.5)*π)    # π the constant!
    if n>2
        π_n   = 1 - ω - r * d
        if n==4
            log_p = log.(p)
        end
    end
    
    if n==2
        log_ω, tan_λ
    elseif n==3
        log_ω, tan_λ, π_n
    else
        log_ω, tan_λ, π_n, log_p
    end
end

convert (generic function with 1 method)

restore

function restore(transformed,r)
    n = length(transformed)
    if n==2
        log_ω, tan_λ = transformed
    elseif n==3
        log_ω, tan_λ, π_n = transformed
    else
        log_ω, tan_λ, π_n, log_p = transformed
    end
    ω = exp.(log_ω)
    λ = atan.(tan_λ)/π + 0.5   # π the constant!
    if n>2
        d = (1 - ω - π_n)/r
        if n==4
            p = atan.(log_p) 
        end
    end
    if n==2
        ω, λ
    elseif n==3
        ω, λ, d
    else
        ω, λ, d, p
    end
end

restore (generic function with 1 method)

Define ODEProblem and solve

# Sample paths of the Keen model
π_bar_n_K = κ_inv(ν * (α + β + δ));
d_bar_K = (κ(π_bar_n_K) - π_bar_n_K)/(α + β);
ω_bar_K = 1 - π_bar_n_K - r * d_bar_K;
λ_bar_K = Φ_inv(α);

Initial values

ω0 = 0.75
λ0 = 0.75
d0 = 0.1
Y0 = 100
T = 300;

ODEProblem

u0 = [-0.288, 1.0, 0.247]
tspan = (0.0,300.0)
p = [α, β, δ, ν, r, (ϕ0,ϕ1), (κ0,κ1,κ2)]
prob = ODEProblem(keen!,u0,tspan,p)

ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 300.0)
u0: [-0.288, 1.0, 0.247]

solve

sol = solve(prob,reltol=1e-6);

Convert to DataFrame

https://diffeq.sciml.ai/stable/features/io/#io

df = DataFrame(sol) # It's the wrong way around!!!!!!!

3 rows × 703 columns (omitted printing of 695 columns)

	x1	x2	x3	x4	x5	x6	x7	x8
	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	-0.288	-0.293876	-0.305295	-0.320146	-0.330904	-0.33368	-0.322919	-0.312928
2	1.0	1.12342	1.46596	2.41279	4.98538	11.5144	35.2791	49.8736
3	0.247	0.24955	0.254226	0.259767	0.263179	0.262936	0.254095	0.246677

State variables

• $ω$ - wage share of the economy
• $λ$ - employment rate
• $d$ - debt (or reciprocal debt: $u=1/d$)
• $Y$ - total output

#state_variables = map(u -> restore(u,r),sol.u);

state_variables = [(restore(u,r),t) for (u,t) in zip(sol.u,sol.t)];

state_variables[2]

((0.7461811191854018, 0.7649302031586827, 0.15779923733286), 0.0748259153505046)

plot([u for (u,t) in state_variables])
xaxis!("ω - wage share of the economy")
yaxis!("λ - employment rate")
#zaxis!("d - debt")

Experiments munging the result

[[u1,u2,u3] for ((u1,u2,u3),t) in state_variables]

1330-element Array{Array{Float64,1},1}:
 [0.7497615922390413, 0.75, 0.10794692536529131]
 [0.7461811191854018, 0.7649302031586827, 0.15779923733286]
 [0.7398108038589221, 0.7945265193981843, 0.25112337052938993]
 [0.7321778575222391, 0.8358964387697227, 0.37057465404153334]
 [0.7263666965926021, 0.8727191304039716, 0.46742844682300577]
 [0.7211036452946267, 0.9117273336598325, 0.5615264153944218]
 [0.7182526371418824, 0.937183540173052, 0.618696007939749]
 [0.7166849652319078, 0.9558377532286569, 0.6586181390481579]
 [0.7162133455271785, 0.9690755840265715, 0.6859638726882704]
 [0.716804590918236, 0.978543607812272, 0.7050051718226277]
 [0.7187297827638306, 0.9853306659147664, 0.7183537329321847]
 [0.7229203640146292, 0.9902674385001823, 0.7278498182572812]
 [0.7317752387441497, 0.9937218650555448, 0.734298880933169]
 ⋮
 [0.8356611821233333, 0.9685433666896768, 0.07177778586506256]
 [0.8357215189457186, 0.9687394282011299, 0.07215217655460368]
 [0.8360608350166773, 0.9688064835542749, 0.07224444259273748]
 [0.8363380269929497, 0.9686563527743229, 0.07187609281675093]
 [0.8362469897749424, 0.9684765059285643, 0.07146073353466269]
 [0.8359183191543937, 0.9684329145945476, 0.07133940355770937]
 [0.8356860288883594, 0.9686164623565069, 0.07169023015628033]
 [0.8358384431155644, 0.9687644979956029, 0.07196671370575947]
 [0.8361799095838119, 0.9687574111957825, 0.07190500575027514]
 [0.8363047698281344, 0.9685635581534429, 0.07145091259050755]
 [0.836100270297997, 0.968449910332632, 0.07118452254630536]
 [0.836006694392986, 0.9684444145580585, 0.07116563305258432]

sol.t

1330-element Array{Float64,1}:
   0.0
   0.0748259153505046
   0.20933786085606826
   0.3735054100661737
   0.5020096204826823
   0.6244492914795371
   0.6984427361334352
   0.7504730619065657
   0.7867571720102611
   0.8128332086393241
   0.8321016278375952
   0.8471244067318621
   0.8592935855891801
   ⋮
 294.19219816611565
 294.71387668708013
 295.25699018406544
 295.84758602400206
 296.37114968486003
 296.92008543869275
 297.58685767529755
 298.0933298151792
 298.67555456940664
 299.32124325626256
 299.83609371027626
 300.0

plot(sol.t, [[u1,u2,u3] for ((u1,u2,u3),t) in state_variables])

BoundsError: attempt to access 3-element Array{Float64,1} at index [1:1330]

Stacktrace:
 [1] throw_boundserror(::Array{Float64,1}, ::Tuple{UnitRange{Int64}}) at .\abstractarray.jl:537
 [2] checkbounds at .\abstractarray.jl:502 [inlined]
 [3] getindex(::Array{Float64,1}, ::UnitRange{Int64}) at .\array.jl:794
 [4] plotly_series_segments(::Plots.Series, ::Dict{Symbol,Any}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Tuple{Float64,Float64}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:724
 [5] plotly_series(::Plots.Plot{Plots.PlotlyBackend}, ::Plots.Series) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:545
 [6] plotly_series(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:862
 [7] plotly_series_json at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:868 [inlined]
 [8] plotly_html_body(::Plots.Plot{Plots.PlotlyBackend}, ::Nothing) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:908
 [9] plotly_html_body at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:903 [inlined]
 [10] html_body at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:873 [inlined]
 [11] standalone_html(::Plots.Plot{Plots.PlotlyBackend}; title::String) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\web.jl:7
 [12] standalone_html(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\web.jl:7
 [13] _show at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:944 [inlined]
 [14] show(::Base.GenericIOBuffer{Array{UInt8,1}}, ::MIME{Symbol("text/html")}, ::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\output.jl:215
 [15] sprint(::Function, ::MIME{Symbol("text/html")}, ::Vararg{Any,N} where N; context::Nothing, sizehint::Int64) at .\strings\io.jl:105
 [16] sprint at .\strings\io.jl:101 [inlined]
 [17] _ijulia_display_dict(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\ijulia.jl:56
 [18] display_dict(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\init.jl:73
 [19] #invokelatest#1 at .\essentials.jl:712 [inlined]
 [20] invokelatest at .\essentials.jl:711 [inlined]
 [21] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\ibuckley\.julia\packages\IJulia\e8kqU\src\execute_request.jl:112
 [22] #invokelatest#1 at .\essentials.jl:712 [inlined]
 [23] invokelatest at .\essentials.jl:711 [inlined]
 [24] eventloop(::ZMQ.Socket) at C:\Users\ibuckley\.julia\packages\IJulia\e8kqU\src\eventloop.jl:8
 [25] (::IJulia.var"#15#18")() at .\task.jl:358

Calculate $Y$

Y_output = Y0*yK(:,2)/𝜆0.*exp((𝛼+𝛽)*tK);

Error using eval
Undefined function 'convert' for input arguments of type 'double'.

Figure 3

# Figure 3
#figure(1)
#plot(yK(:,1), yK(:,2)), hold on;
#plot(bar_omega_K, bar_lambda_K, 'ro')

Figure 4

#% Figure 4
#figure(2)
#yyaxis right
#plot(tK, yK(:,1),tK,yK(:,2),tK,yK(:,3))
#ylabel('\omega,\lambda,d')

#yyaxis left
#plot(tK, Y_output)
#ylabel('Y')

#title(['\omega_0 = ', num2str(omega0,txt_format), ...
#    ', \lambda_0 = ', num2str(lambda0, txt_format), ...
#    ', d_0 = ', num2str(d0, txt_format), ...
#    ', Y_0 = ', num2str(Y0, txt_format)])
#legend('Y','\omega', '\lambda','d')
#print('keen_time_convergent','-depsc')


#omega0 = 0.75;
#lambda0 = 0.70;
#d0 = 0.1;
#Y0 = 100;
#T = 300;

#[tK,yK] = ode15s(@(t,y) keen(y), [0 T], convert([omega0, lambda0, d0]), options);
#yK = retrieve(yK);
#Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);

Figure 5

#% Figure 5
#figure(3)
#yyaxis right
#plot(tK, yK(:,1),tK,yK(:,2))
#ylabel('\omega,\lambda')

#yyaxis left
#plot(tK,yK(:,3),tK, Y_output)
#ylabel('d,Y')

#title(['\omega_0 = ', num2str(omega0,txt_format), ...
#    ', \lambda_0 = ', num2str(lambda0, txt_format), ...
#    ', d_0 = ', num2str(d0, txt_format), ...
#    ', Y_0 = ', num2str(Y0, txt_format)])
#legend('d', 'Y','\omega','\lambda')
#print('keen_time_divergent','-depsc')

Auxiliary functions

Do not run these in the notebook. MATLAB does not support functions outside of files of the same name.

### Auxiliary functions
## Instead of integrating the system in terms of omega, lambda and d, we
## decided to use:
##
## log_omega = log(omega), 
## tan_lambda = tan((lambda-0.5)*pi), and 
## pi_n=1-omega-r*d
## log_p = log(p)
##
## Experiments show that the numerical methods work better with
## these variables.

function new = convert(old)
    # This function converts
    # [omega, lambda, d, p]
    # to
    # [log_omega, tan_lambda, pi_n, log_p]
    # It also work with the first 2 or 3 elements alone.
    # each of these variables might also be a time-dependent vector
    
    n = size(old,2); #number of variables
    new = zeros(size(old));
    
    new(:,1) = log(old(:,1));
    new(:,2) = tan((old(:,2)-0.5)*pi);
    if n>2
        new(:,3) = 1-old(:,1)-r*old(:,3);
        if n==4
            new(:,4) = log(old(:,4));    
        end
    end
end

Error: Function definition not supported in this context. Create functions in code file.

function old = retrieve(new)
    # This function retrieves [omega, lambda, d, p]
    # from
    # [log_omega, tan_lambda, pi_n, log_p]
    #  
    #  It also work with the first 2 or 3 elements alone.
    # each of these variables might also be a time-dependent vector
    
    n = size(new,2); #number of variables    
    old = zeros(size(new));
    
    old(:,1) = exp(new(:,1));
    old(:,2) = atan(new(:,2))/pi+0.5;
    if n>2
        old(:,3) = (1-old(:,1)-new(:,3))/r;
        if n==4
            old(:,4) = exp(new(:,4));
        end
    end
end

Error: Function definition not supported in this context. Create functions in code file.

function f = keen(y)    
    f = zeros(3,1);
    log_omega = y(1);
    tan_lambda = y(2);
    pi_n = y(3);
    lambda = atan(tan_lambda)/pi+0.5;
    omega = exp(log_omega);
    g_Y = fun_kappa(pi_n)/nu-delta;    
    
    % Equation (47)
    f(1) = fun_Phi(lambda)-alpha;                                   %d(log_omega)/dt
    f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta);             %d(tan_lambda)/dt
    f(3) = -omega*f(1)-r*(fun_kappa(pi_n)-pi_n)+(1-omega-pi_n)*g_Y; %d(pi_n)/dt
end

Error: Function definition not supported in this context. Create functions in code file.

Test plotting

https://github.com/imatlab/imatlab
https://github.com/plotly/plotly-graphing-library-for-matlab

tx = linspace (-8, 8, 41);
ty = tx;
[xx, yy] = meshgrid (tx, ty);
r = sqrt (xx .^ 2 + yy .^ 2) + eps;
tz = sin (r) ./ r;

mesh(tx, ty, tz)

Drawing...